Purpose: To measure the capacity for each strain to grow as a biofilm
Protocol:
library(tidyverse)
library(xlsx)
library(broom)
library(rstatix)
library(ape)
library(ggbeeswarm)
library(formattable)
library(ggpubr)
library(cowplot)
library(corrplot)
library(ggpubr)
library(grid)
library(gridGraphics)
`%!in%` <- negate(`%in%`)
Short cuts:
strains <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35", "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2", "UM224", "Gv5-1", "C0179B4", "C0056B5", "Gv101", "C0100B2", "C0101A1", "CMW7778B")
strains0 <- c("14018", "315-A", "C0011E4", "JCP7275", "C0084H9", "JCP8108", "JCP8017B", "C0040C2" , "UM35", "JCP8066", "C0093B3", "AMD", "C0096A1", "C0102A2", "UM224", "5-1", "C0179B4", "C0056B5", "101", "C0100B2", "C0101A1", "CMW7778B")
species <- c("Gardnerella vaginalis", "Gardnerella sp. 2", "Gardnerella sp. 3", "Gardnerella piotii", "Gardnerella leopoldii", "Gardnerella swidsinskii", "Gardnerella sp. 7", "Gardnerella sp. 8", "Gardnerella sp. 10", "Gardnerella sp. 11", "Gardnerella sp. 12")
figureOut <- "../../experiments_figures"
suppTableOut <- "../../experiments_figures/supplementary_posthoc_tables.xlsx"
Import data
midlog_raw_data <- read_csv("./data/biofilm_midlog_raw_data.csv")
rawPlanktonicODs <- read_csv("./data/planktonic_od_raw_data.csv")
rawBiofilmODs <- read_csv("./data/biofilm_od_raw_data.csv")
rawSafPlanktonicODs <- read_csv("./data/safranin_planktonic_raw_data.csv")
rawSafBiofilmODs <- read_csv("./data/safranin_biofilm_raw_data.csv")
Plot settings
theme_set(theme_bw()+
theme(panel.grid = element_blank()))
speciesColor <- scale_color_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73",
"darkgreen", "gold3", "darkred", "darkblue",
"mediumpurple4"))
speciesFill <- scale_fill_manual(values=c("#CC79A7", "#D55E00", "#56B4E9", "#0072B2", "#E69F00", "#009E73",
"darkgreen", "gold3", "darkred", "darkblue",
"mediumpurple4"))
Propagate twice and then incubate to mid log. 1:10 Dilution for mid Log: Protocol 1: 4 hour incubation Protocol 2: 7 hour incubation Protocol 3: 15 hour incubation ## Clean data
setupDF <- midlog_raw_data %>%
group_by(Number, BioRep, Species, Strain, Batch, Step) %>%
mutate(OD=OD-mean(.$OD[.$Number=="Blank"]),
OD=case_when(OD>0~OD,
OD<=0~0)) %>%
filter(Number!="Blank") %>%
mutate(Step=factor(Step, levels=c("Overnight 1", "Overnight 2", "Mid Log"))) %>%
mutate(Strain=factor(Strain, levels=strains),
Species=factor(Species, levels=species),
Sample=paste(Strain, Batch, BioRep, sep="."),
BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
BatchRep=paste(BatchNum, BioRep, sep="."),
StrainN=as.numeric(Strain),
rectFill=case_when(StrainN %% 2 == 0 ~ "A",
StrainN %% 2 != 0 ~ "B")) %>%
ungroup %>%
with_groups(c("Sample", "Number", "BioRep", "Batch", "BatchNum", "BatchRep", "StrainN", "rectFill", "Species", "Strain", "Protocol", "Step"), summarize, ODm=mean(OD), ODsd=sd(OD))
Assess mid log growth for removing poor growers with OD600 < 0.05
setupDF %>%
filter(Step=="Mid Log") %>%
ggplot(aes(x=Strain, y=ODm)) +
geom_point(aes(shape=BioRep)) +
geom_hline(aes(yintercept=0.05), linetype=2, color="gray") +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=bquote(OD[600]), shape="Biological Replicate")
Create quality filter for OD600 > 0.05
midLogFail <- setupDF %>%
filter(Step=="Mid Log",
ODm < 0.05) %>%
.$Sample
midLogFail
## [1] "14018.1A1.A" "14018.1A1.B" "14018.1A2.A" "14018.1A2.B"
## [5] "AMD.2A1.B" "C0096A1.2A2.A" "C0096A1.2A3.A" "C0100B2.1AB3.A"
## [9] "C0100B2.1AB3.B" "C0102A2.1A1.B" "C0102A2.1A2.B" "CMW7778B.2A1.A"
## [13] "CMW7778B.2A1.B" "CMW7778B.2A2.A" "CMW7778B.2A2.B" "CMW7778B.2A3.A"
## [17] "Gv5-1.1AB3.A" "Gv5-1.1AB3.B" "Gv5-1.1B1.A" "Gv5-1.1B1.B"
## [21] "Gv5-1.1B2.B" "JCP7275.1AB3.A" "JCP7275.1AB3.B" "JCP7275.1AB4.B"
## [25] "JCP7275.1B1.A" "JCP7275.1B1.B" "JCP7275.1B2.A" "JCP7275.1B2.B"
## [29] "JCP8066.1A2.A" "JCP8066.1A2.B" "JCP8066.1AB3.A" "JCP8066.1AB3.B"
## [33] "JCP8108.1B1.A" "JCP8108.1B2.A" "JCP8108.1B2.B" "UM224.1A1.A"
## [37] "UM224.1A2.B" "UM35.3A1.B" "UM35.3A2.A" "UM35.3A2.B"
## [41] "UM35.3A3.B"
overnightPlot <- setupDF %>%
filter(Sample %!in% midLogFail) %>%
ggplot(aes(x=Strain, y=ODm)) +
geom_point(aes(shape=BatchNum), alpha=0) +
geom_rect(aes(xmin=StrainN-.5, xmax=StrainN+.5, ymin = 0, ymax = 1, fill=rectFill), alpha = 0.2, show.legend = FALSE) +
geom_pointrange(aes(ymin=(ODm-ODsd), ymax=(ODm+ODsd), shape=BatchNum, color=Species), position = position_quasirandom(), size=0.3, alpha=0.7, show.legend = FALSE) +
speciesColor +
scale_fill_manual(values=c("white", "gray"), labels=c("A","B")) +
ylim(0,1) +
facet_grid(.~Step) +
coord_flip() +
scale_x_discrete(limits=rev) +
theme(axis.text.x = element_text(size=8),
axis.text.y = element_text(size=8),
axis.title = element_text(size=11),
panel.grid = element_blank(),
legend.position = "none") +
labs(x=NULL, y=bquote(OD[600]), shape="Experiment Batch")
#annotation bars
protocolBar <- setupDF %>%
select(Strain, Species, Protocol) %>%
mutate(Protocol=as.character(Protocol)) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
scale_x_discrete(limits=rev) +
scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
coord_flip() +
theme_void() +
theme(legend.position = "none")
speciesBar <- setupDF %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
scale_x_discrete(limits=rev) +
speciesFill +
coord_flip() +
theme_void() +
theme(legend.position = "none")
#Legends
protocol_legend <- ggpubr::get_legend(protocolBar +
guides(color=guide_legend(ncol = 1)) +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
species_legend <- ggpubr::get_legend(speciesBar +
guides(color=guide_legend(ncol = 1)) +
theme(legend.position = "bottom",
legend.direction = "horizontal"))
experiment_legend <- ggpubr::get_legend(
overnightPlot +
guides(shape=guide_legend(nrow = 1, override.aes = list(alpha=1, size=4))) +
theme(legend.position = "bottom"))
a <- plot_grid(overnightPlot, speciesBar, protocolBar, nrow = 1, align = "h", axis = "bt", rel_widths = c(1, 0.02, 0.02))
b <- plot_grid(a, as_ggplot(species_legend), nrow=2, rel_heights = c(1, 0.3))
c <- plot_grid(as_ggplot(protocol_legend), as_ggplot(experiment_legend), nrow=1)
plot_grid(b, c, ncol=1, nrow=2, rel_heights = c(1, 0.1))
planktonicODs <- rawPlanktonicODs %>%
group_by(Batch) %>%
mutate(planktonicOD=OD-mean(.$OD[.$Number=="Blank"]),
planktonicOD=case_when(planktonicOD>0~planktonicOD,
planktonicOD<=0~0)) %>%
filter(Number!="Blank") %>%
select(-OD, -Measure) %>%
ungroup
biofilmODs <- rawBiofilmODs %>%
group_by(Batch) %>%
mutate(biofilmOD=OD-mean(.$OD[.$Number=="Blank"]),
biofilmOD=case_when(biofilmOD>0~biofilmOD,
biofilmOD<=0~0)) %>%
filter(Number!="Blank") %>%
select(-OD, -Measure) %>%
ungroup
percentOD <- planktonicODs %>%
left_join(biofilmODs) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
filter(Sample %!in% midLogFail) %>% # quality filter for poor growth at mid log phase
mutate(edgeWell=str_detect(well, pattern="[A|H]"),
edgeWell=factor(edgeWell, levels=c(TRUE, FALSE)),
innerWell=str_detect(well, pattern="[D|E]"),
innerWell=factor(innerWell, levels=c(TRUE, FALSE)),
topWell=str_detect(well, pattern="A"),
topWell=factor(topWell, levels=c(TRUE, FALSE)))
head(percentOD)
## # A tibble: 6 × 13
## well Number BioRep Strain Species Protocol Batch planktonicOD biofilmOD
## <chr> <chr> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl>
## 1 A2 1 A C0011E4 Gardnerella… 1 1A1 0.196 0.0663
## 2 B2 1 A C0011E4 Gardnerella… 1 1A1 0.193 0.0773
## 3 C2 1 A C0011E4 Gardnerella… 1 1A1 0.172 0.0833
## 4 D2 1 A C0011E4 Gardnerella… 1 1A1 0.138 0.0913
## 5 E2 1 B C0011E4 Gardnerella… 1 1A1 0.0895 0.152
## 6 F2 1 B C0011E4 Gardnerella… 1 1A1 0.0845 0.221
## # ℹ 4 more variables: Sample <chr>, edgeWell <fct>, innerWell <fct>,
## # topWell <fct>
edgeWellPlot <- percentOD %>%
mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
Strain=factor(Strain, levels = strains),
Species=factor(Species, levels = species)) %>%
ggplot(aes(x=Strain, y=percentBiofilm, shape=BioRep, color=edgeWell)) +
geom_quasirandom() +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "bottom",
legend.direction = "vertical") +
labs(x=NULL, y="Percent Biofilm", shape="Biological Replicate", title="Edge Wells")
innerWellPlot <- percentOD %>%
mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
Strain=factor(Strain, levels = strains),
Species=factor(Species, levels = species)) %>%
ggplot(aes(x=Strain, y=percentBiofilm, shape=BioRep, color=innerWell)) +
geom_quasirandom() +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "bottom",
legend.direction = "vertical") +
labs(x=NULL, y=NULL, shape="Biological Replicate", title="Inner Wells")
plot_grid(edgeWellPlot, innerWellPlot, nrow=1)
percentODDF0 <- percentOD %>%
mutate(percentBiofilm=(biofilmOD/(biofilmOD+planktonicOD))*100,
Strain=factor(Strain, levels = strains),
Species=factor(Species, levels = species),
Protocol=as.character(Protocol),
BatchRep=paste(Batch, BioRep, sep=".")) %>%
gather("measure", "value", c(planktonicOD, biofilmOD, percentBiofilm)) %>%
mutate(measure=factor(measure, levels=c("biofilmOD", "planktonicOD", "percentBiofilm"), labels=c("Biofilm OD", "Planktonic OD", "Percent OD")),
BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
BatchRep=paste(BatchNum, BioRep, sep="."))
percentODDF0 %>%
ggplot(aes(x=Strain, y=value)) +
geom_boxplot(alpha=0) +
geom_quasirandom(aes(color=edgeWell), alpha=0.7) +
scale_y_continuous(position = "right") +
facet_grid(measure~., scales="free") +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=NULL, shape="Biological Replicate")
Remove Edge Wells
naSamples <- percentODDF0 %>%
filter(edgeWell==FALSE) %>%
spread(measure, value) %>%
filter(is.na(`Percent OD`)) %>%
.$Sample %>%
unique
naSamples
## [1] "C0093B3.3A2.A" "C0093B3.3A2.B"
percentODDF <- percentODDF0 %>%
filter(edgeWell==FALSE, # remove samples
Sample %!in% naSamples) %>% # remove samples with no growth
mutate(method="Percent OD")
non4StrainsPOD <- percentODDF0 %>%
mutate(Sample=paste(Number, Sample, sep=".")) %>%
with_groups(Strain, summarize, n=n_distinct(Sample)) %>%
filter(n!=4) %>%
.$Strain
percentODDF %>%
filter(Strain %in% non4StrainsPOD) %>%
mutate(Sample=paste(Number, Sample, sep=".")) %>%
with_groups(Strain, mutate, n=n_distinct(Sample)) %>%
with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "measure", "BatchNum", "n"), summarize, valueM=mean(value), valueSD=sd(value)) %>%
mutate(labelYpos=case_when(measure=="Biofilm OD"~0.75,
measure=="Planktonic OD"~0.75,
measure=="Percent OD"~110)) %>%
ggplot() +
geom_boxplot(aes(x=Strain, y=valueM), alpha=0) +
geom_pointrange(aes(x=Strain, y=valueM, ymin=valueM-valueSD, ymax=valueM+valueSD, shape=BatchRep), alpha=0.7, size=0.25, position=position_quasirandom()) +
geom_text(aes(x=Strain, y=labelYpos, label=n)) +
scale_y_continuous(position = "right") +
facet_grid(measure~., scales="free") +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=NULL)
## annotation portions of plot
#protocol bar
protocolBar <- percentODDF %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=2, fill=Protocol), stat="identity", width = 1.1) +
scale_fill_manual(values=c("cornflowerblue", "burlywood1", "forestgreen")) +
theme_void() +
theme(legend.position = "none")
# species bar
speciesBar <- percentODDF %>%
select(Strain, Species, Protocol) %>%
unique %>%
ggplot() +
geom_bar(aes(x=Strain, y=1, fill=Species), stat="identity", width=1.1) +
speciesFill +
theme_void() +
theme(legend.position = "none")
percentODPlot <- percentODDF %>%
with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "measure", "BatchNum"), summarize, valueM=mean(value), valueSD=sd(value)) %>%
ggplot() +
geom_boxplot(aes(x=Strain, y=valueM), alpha=0) +
geom_pointrange(aes(x=Strain, y=valueM, ymin=valueM-valueSD, ymax=valueM+valueSD, color=Species), alpha=0.7, size=0.25, position=position_quasirandom()) +
speciesColor +
facet_grid(measure~., scales="free", switch = "y") +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none",
strip.placement = "outside") +
labs(x=NULL, y=NULL)
ODPlots <- plot_grid(protocolBar, speciesBar, percentODPlot, ncol = 1, rel_heights = c(0.06, 0.06, 1), align = "v", axis = "lr")
ODPlots
strain_OD_DF <- percentODDF %>%
spread(measure, value) %>%
select(-`Biofilm OD`, -`Planktonic OD`) %>%
dplyr::rename(percentOD=`Percent OD`) %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mPercentOD=mean(percentOD), sdPercentOD=sd(percentOD))
strain_OD_DF %>%
ggplot(aes(sample=mPercentOD)) +
geom_qq_line() +
geom_qq()
#### One-way ANOVA
strain_OD_DF %>%
anova_test(mPercentOD~Strain) %>%
formattable
| Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|
| Strain | 20 | 55 | 5.49 | 2.6e-07 |
|
0.666 |
ODPostHoc <- strain_OD_DF %>%
tukey_hsd(mPercentOD~Strain, p.adjust.method = "BH")
# ODPostHoc %>%
# arrange(p.adj) %>%
# write.xlsx(file=suppTableOut, sheetName="TableS5_percent_od_posthoc", col.names=TRUE, append=TRUE)
pvals <- ODPostHoc %>%
select(group1, group2, p.adj) %>%
dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix()
ODPostHoc %>%
select(group1, group2, estimate) %>%
dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix() %>%
corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)
grid.echo()
percODPostHocPlot <- grid.grab()
Percent biofilm of Strains
strain_OD_DF %>%
with_groups(Strain, summarize, `Percent OD Biofilm`=mean(mPercentOD)) %>%
arrange(-`Percent OD Biofilm`) %>%
formattable()
| Strain | Percent OD Biofilm |
|---|---|
| C0040C2 | 97.58907 |
| JCP8066 | 91.16070 |
| C0056B5 | 90.17673 |
| C0096A1 | 84.37130 |
| C0100B2 | 84.26520 |
| C0102A2 | 83.82580 |
| C0179B4 | 81.63978 |
| Gv5-1 | 80.78987 |
| UM224 | 80.75329 |
| AMD | 76.27148 |
| UM35 | 70.98144 |
| C0101A1 | 69.57704 |
| JCP8017B | 64.14116 |
| 315-A | 61.51312 |
| JCP8108 | 59.02140 |
| C0011E4 | 51.01499 |
| C0084H9 | 50.11188 |
| Gv101 | 48.29067 |
| C0093B3 | 41.46202 |
| JCP7275 | 36.34031 |
| CMW7778B | 27.52904 |
Average of all strains
strain_OD_DF %>%
summarize(`Average Percent OD Biofilm`=mean(mPercentOD)) %>%
formattable()
| Average Percent OD Biofilm |
|---|
| 70.12149 |
Load tree
We choose the weights as \(_{wij} = 1/d_{ij}\) , where the d’s is the distances measured on the tree:
We can now perform the analysis with Moran’s I:
safPlankODDF <- rawSafPlanktonicODs %>%
group_by(Batch) %>%
mutate(safPlanktonicOD=OD-mean(.$OD[.$Number=="Blank"]),
OD=case_when(OD>0~OD,
OD<=0~0)) %>%
filter(Number!="Blank") %>%
ungroup %>%
mutate(row=str_extract(well, pattern="[A-H]")) %>%
select(-OD, -Measure, -well)
safSafODDF <- rawSafBiofilmODs %>%
group_by(Batch) %>%
mutate(safraninOD=OD-mean(.$OD[.$Number=="Blank"]),
OD=case_when(OD>0~OD,
OD<=0~0)) %>%
filter(Number!="Blank") %>%
ungroup %>%
mutate(row=str_extract(well, pattern="[A-H]")) %>%
select(-OD, -Measure, -well)
safraninOD1 <- safPlankODDF %>%
full_join(safSafODDF) %>%
select(row, everything(), safPlanktonicOD, safraninOD) %>%
mutate(Sample=paste(Strain, Batch, BioRep, sep=".")) %>%
filter(Sample %!in% midLogFail) %>% # quality filter for poor growth at mid log phase
mutate(edgeWell=str_detect(row, pattern="[A|H]"),
edgeWell=factor(edgeWell, levels=c(TRUE, FALSE)),
innerWell=str_detect(row, pattern="[D|E]"),
innerWell=factor(innerWell, levels=c(TRUE, FALSE)),
topWell=str_detect(row, pattern="A"),
topWell=factor(topWell, levels=c(TRUE, FALSE)))
head(safraninOD1)
## # A tibble: 6 × 13
## row Number BioRep Strain Species Protocol Batch safPlanktonicOD safraninOD
## <chr> <chr> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl>
## 1 A 1 A C0011E4 Gardner… 1 1A1 0.151 0.202
## 2 B 1 A C0011E4 Gardner… 1 1A1 0.142 0.217
## 3 C 1 A C0011E4 Gardner… 1 1A1 0.136 0.232
## 4 D 1 A C0011E4 Gardner… 1 1A1 0.159 0.201
## 5 E 1 B C0011E4 Gardner… 1 1A1 0.152 0.179
## 6 F 1 B C0011E4 Gardner… 1 1A1 0.116 0.103
## # ℹ 4 more variables: Sample <chr>, edgeWell <fct>, innerWell <fct>,
## # topWell <fct>
Organize data and create dataframe
safDF <- safraninOD1 %>%
mutate(Strain=factor(Strain, levels = strains),
Species=factor(Species, levels = species),
Protocol=as.character(Protocol),
BatchNum=str_extract(Batch, "(?<=[A|B])[0-9]"),
BatchRep=paste(BatchNum, BioRep, sep="."),
method="Safranin Stain") %>%
filter(Sample %!in% naSamples) # remove samples with no growth
Check if edge wells should be removed
safDF %>%
ggplot(aes(x=Strain, y=safraninOD)) +
geom_boxplot(alpha=0) +
geom_quasirandom(aes(color=edgeWell), alpha=0.7) +
scale_y_continuous(position = "right") +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=NULL, shape="Biological Replicate")
Do not remove edge wells ## Extra biological replicates
non4StrainsSaf <- safDF %>%
mutate(Sample=paste(Number, Sample, sep=".")) %>%
with_groups(Strain, summarize, n=n_distinct(Sample)) %>%
filter(n!=4) %>%
.$Strain
safDF %>%
filter(Strain %in% non4StrainsSaf) %>%
mutate(Sample=paste(Number, Sample, sep=".")) %>%
with_groups(Strain, mutate, n=n_distinct(Sample)) %>%
with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "BatchNum", "n"), summarize, safraninM=mean(safraninOD), safraninSD=sd(safraninOD)) %>%
ggplot() +
geom_boxplot(aes(x=Strain, y=safraninM), alpha=0) +
geom_pointrange(aes(x=Strain, y=safraninM, ymin=safraninM-safraninSD, ymax=safraninM+safraninSD, shape=BatchRep), alpha=0.7, position=position_quasirandom(), size=0.3) +
geom_text(aes(x=Strain, y=0.575, label=n)) +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Batch")
safPlot <- safDF %>%
ggplot(aes(x=Strain, y=safraninOD, color=Species, shape=BatchRep)) +
geom_quasirandom(alpha=0.75) +
theme(axis.text.x = element_text(angle=45, hjust=1)) +
labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Biological Replicate")
safPlot2 <- safDF %>%
with_groups(c("Number", "Strain", "Species", "Protocol", "Batch", "method", "Sample", "BatchRep", "BatchNum"), summarize, safraninM=mean(safraninOD), safraninSD=sd(safraninOD)) %>%
ggplot() +
geom_boxplot(aes(x=Strain, y=safraninM), alpha=0) +
geom_pointrange(aes(x=Strain, y=safraninM, ymin=safraninM-safraninSD, ymax=safraninM+safraninSD, color=Species), alpha=0.7, position=position_quasirandom(), size=0.3) +
speciesColor +
theme(axis.text.x = element_text(angle=45, hjust=1),
legend.position = "none") +
labs(x=NULL, y=bquote("Safranin OD"[450]), shape="Batch")
safPlots <- plot_grid(protocolBar, speciesBar, safPlot2, ncol = 1, rel_heights = c(0.06, 0.06, 1), align = "v", axis = "lr")
safPlots
safDF %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSafOD=mean(safraninOD), sdSafOD=sd(safraninOD)) %>%
ggplot(aes(sample=mSafOD)) +
geom_qq_line() +
geom_qq()
#### One-way ANOVA
strain_saf_DF <- safDF %>%
with_groups(c(Sample, Species, Strain, Batch, BioRep), summarize, mSafOD=mean(safraninOD), sdSafOD=sd(safraninOD))
strain_saf_DF %>%
anova_test(mSafOD~Strain) %>%
formattable
| Effect | DFn | DFd | F | p | p<.05 | ges |
|---|---|---|---|---|---|---|
| Strain | 20 | 55 | 4.469 | 5.31e-06 |
|
0.619 |
safPostHoc <- strain_saf_DF %>%
tukey_hsd(mSafOD~Strain, p.adjust.method = "BH")
# safPostHoc %>%
# arrange(p.adj) %>%
# write.xlsx(file=suppTableOut, sheetName="TableS6_safranin_posthoc", col.names=TRUE, append=TRUE)
pvals <- safPostHoc %>%
select(group1, group2, p.adj) %>%
dplyr::rename(var1=group2, var2=group1, cor=p.adj) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix()
safPostHoc %>%
select(group1, group2, estimate) %>%
dplyr::rename(var1=group2, var2=group1, cor=estimate) %>%
cor_spread() %>%
column_to_rownames(var="rowname") %>%
as.matrix() %>%
corrplot(is.corr=FALSE, type = "lower", tl.col="black", p.mat = pvals, insig = "label_sig", sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9)
grid.echo()
safraninPostHocPlot <- grid.grab()
We can now perform the analysis with Moran’s I:
plot_grid(protocolBar, speciesBar, percentODPlot + speciesColor + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()), safPlot2, ncol = 1, rel_heights = c(0.07, 0.07, 1, 0.5), align = "v", axis = "lr", labels = c("A", "", "", "B"), label_size = 15)
#ggsave(file.path(figureOut, paste(Sys.Date(), "Figure4_strain_biofilm.png", sep="_")))
plot_grid(ggplot()+geom_blank(), percODPostHocPlot, ggplot()+geom_blank(), safraninPostHocPlot, nrow = 1, rel_widths = c(0.1, 1, 0.1, 1), labels = c(" ", "A", " ", "B"), label_size = 20)
#ggsave(file.path(figureOut, paste(Sys.Date(), "FigureS3_biofilm_posthocs.png", sep = "_")))
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gridGraphics_0.5-1 corrplot_0.92 cowplot_1.1.3 ggpubr_0.6.0
## [5] formattable_0.2.1 ggbeeswarm_0.7.2 ape_5.8 rstatix_0.7.2
## [9] broom_1.0.6 xlsx_0.6.5 lubridate_1.9.3 forcats_1.0.0
## [13] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [17] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.5 beeswarm_0.4.0 xfun_0.46 bslib_0.7.0
## [5] htmlwidgets_1.6.4 rJava_1.0-11 lattice_0.22-6 tzdb_0.4.0
## [9] vctrs_0.6.5 tools_4.4.0 generics_0.1.3 parallel_4.4.0
## [13] fansi_1.0.6 highr_0.11 pkgconfig_2.0.3 lifecycle_1.0.4
## [17] farver_2.1.2 compiler_4.4.0 munsell_0.5.1 carData_3.0-5
## [21] vipor_0.4.7 htmltools_0.5.8.1 sass_0.4.9 yaml_2.3.9
## [25] crayon_1.5.3 pillar_1.9.0 car_3.1-2 jquerylib_0.1.4
## [29] cachem_1.1.0 abind_1.4-5 nlme_3.1-165 tidyselect_1.2.1
## [33] digest_0.6.36 stringi_1.8.4 labeling_0.4.3 fastmap_1.2.0
## [37] colorspace_2.1-0 cli_3.6.3 magrittr_2.0.3 xlsxjars_0.6.1
## [41] utf8_1.2.4 withr_3.0.0 scales_1.3.0 backports_1.5.0
## [45] bit64_4.0.5 timechange_0.3.0 rmarkdown_2.27 bit_4.0.5
## [49] ggsignif_0.6.4 hms_1.1.3 evaluate_0.24.0 knitr_1.48
## [53] rlang_1.1.4 Rcpp_1.0.13 glue_1.7.0 vroom_1.6.5
## [57] rstudioapi_0.16.0 jsonlite_1.8.8 R6_2.5.1